Crystallization of Carbon Oxygen Mixtures in White Dwarf Stars 
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We determine the piiase diagram for dense carbon/ oxygen mixtures in White Dwarf (WD) star 
interiors using molecular dynamics simulations involving hquid and soHd phases. Our phase diagram 
agrees well with predictions from Ogata et al. and Medin and Gumming and gives lower melting 
temperatures than Segretain et al. Observations of WD crystallization in the globular cluster NGC 
6397 by Winget et al. suggest that the melting temperature of WD cores is close to that for pure 
carbon. If this is true, our phase diagram implies that the central oxygen abundance in these stars is 
less than about 60%. This constraint, along with assumptions about convection in stellar evolution 
models, limits the effective S factor for the ^^0(0,7)^*^0 reaction to Saoo < 170 keV barns. 
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Observations of cooling White Dwarf (WD) stars pro- 
vide important information on the ages of stellar systems 
[T] . The interior of a WD is a coulomb plasma of ions and 
a degenerate electron gas. As the star cools this plasma 
crystallizes. Winget et al. recently observed effects from 
the latent heat of crystallization on the luminosity func- 
tion of WDs in the globular cluster Ngc 6397 [2; . Winget 
et al.'s observations constrain the melting temperature of 
the carbon and oxygen mixtures expected in these WD 
cores. This temperature depends on the ratio of carbon 
to oxygen. Therefore observations of crystallization may 
provide information on WD composition. 

The ratio of carbon to oxygen in WD stars is very in- 
teresting. It depends on the reaction ^^C(a,7)^^0. De- 
spite a great deal of effort, see for example [3], the stellar 
rate for this reaction remains one of the most important 
unsettled rates left in Nuclear Astrophysics ^. Further- 
more, the ratio of carbon to oxygen in massive stars is 
important for their subsequent evolution and nucleosyn- 
thesis [3. Therefore, a measurement of the carbon to 
oxygen ratio in a WD could be very important. 

To determine the C/0 ratio from observations of the 
melting temperature one needs the phase diagram for 
carbon and oxygen mixtures. Segretain et al. calculated 
the phase diagram assuming a local density model for 
the free energy of the sohd ^ . While, Ogata et al. [7] , [S] 
and DeWitt et al. [5], [10] calculated the phase diagram 
based on Monte Carlo or Molecular Dynamics (MD) sim- 
ulation free energies for both the liquid and solid phases. 
Recently Potekhin et al. have made accurate calculations 
of the free energy of liquid mixtures [TT] , [T^ and Medin 
and Gumming calculated the phase diagram for both bi- 
nary mixtures such as C/0 and much more complicated 
multicomponent mixtures |13j . All of these calculations 
are very sensitive to small errors in the free energy differ- 
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ence between liquid and solid phases. Indeed Segretain 
et al. predict higher melting temperatures and a spindle 
type phase diagram while both Ogata et al. and Medin 
and Cumming predict lower melting temperatures and 
an azeotrope type phase diagram. 

In this paper, we present direct two phase molecular 
dynamics simulations of the carbon / oxygen phase dia- 
gram to address these uncertainties. The systematic er- 
rors of our simulations may be different from previous free 
energy calculations. We discuss our formalism, present 
results for the C/0 phase diagram, and present possi- 
ble limits on the central oxygen concentration of WDs in 
NGC 6397 and the effective astrophysical S factor that 
describes the ^^C(a,7)^^0 cross section. 

We describe our two-phase MD simulation formalism. 
This is very similar to what we used earlier for the freez- 
ing of rapid proton capture nucleosynthesis ash on ac- 
creting neutron stars [M]. We assume the electrons 
form a degenerate Fermi gas. The ions are fully pres- 
sure ionized and interact with each other via screened 
Coulomb interactions. The potential between the zth 
and jth ion is assumed to be Vij{r) — ZiZje^e~^^^/r. 
Here the ion charges are Zi and Zj, r is their separation 
and the electron screening length is A. For cold rela- 
tivistic electrons, the Thomas Fermi screening length is 
A""'^ = 2Q^^'^kp j-K^I"^ where the electron Fermi momen- 
tum kp is kp — (37r^ne)^'^ and a is the fine structure 
constant. Finally the electron density n^ is equal to the 
ion charge density, n^ = {Z)n, where n is the ion den- 
sity and {Z) is the average charge. Our simulations are 
classical and for historical reasons we have neglected the 
electron mass. Including the electron mass at a White 
Dwarf central density of a few 10^ g/cm^ will reduce the 
screening length by about 20%. We expect this to have 
only a very small effect on our computed phase diagram, 
however see [TS]. Likewise quantum effects, at these den- 
sities, should also have very small effects on the phase 
diagram because the parameter r<, describing the ratio 
of the ion sphere radius to the ion Bohr radius is large 



rs ~ 18000 [IS], see also [17]. 

We now describe the initial conditions for our classical 
MD simulations. It can be difficult to obtain an equilib- 
rium crystal configuration for a large system involving a 
mixture of ions. Therefore, we start with a very small 
system of 432 ions with random coordinates at a high 
temperature and cool the system a number of times by 
re-scaling the velocities until the system solidifies. Next 
four copies of this solid configuration were placed in the 
top half of a larger simulation volume along with four 
copies of a 432 ion liquid configuration. The resulting 
system with 3456 ions was evolved in time until it fully 
crystalized. Finally, four copies of this 3456 ion crystal 
were placed in the top half of the final simulation volume 
and four copies of a 3456 ion liquid configuration were 
placed in the bottom half. This final system has 27648 
ions and consists of a solid phase above a liquid phase. 
There are two liquid-solid interfaces. The first is near the 
middle of the simulation volume and the second is at the 
top. This is because of the periodic boundary conditions 
that identify the top face with the bottom face. 

The simulations can be characterized by an average 
Coulomb parameter F, 
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Here {Z'°l'^) is an average over the ion charges, T is the 

temperature, and the electron sphere radius a^ is Og = 

(3/47rn,)i/3. 

All of our simulations are run for the same electron 
density of rig — 5.026 x 10"'* fm"-^. Since the pressure 
is dominated by the electronic contribution, constant 
electron density corresponds, approximately, to constant 
pressure. Ignoring quantum effects, the density can be 
scaled to other values by also changing the temperature 
T so that the value of F, see Eq. [l] remains the same. 

We have performed six simulations with parameters as 
indicated in Table llj We evolve the system in time using 
the simple velocity Verlet algorithm |TS] with a time step 
Ai = 25 fm/c for the pure carbon simulation and 100 
fm/c for the five carbon/ oxygen mixture simulations. 
We use periodic boundary conditions. Our simulation 
volume is large enough so that the box length L is much 
larger than the electron screening length A. The ratio 
of the force on two ions separated by a distance L/2, 
compared to the force on two ions separated by the ion 
sphere radius a = (3/47™)!/^ is F{L/2)/F{a) « lO"'^. 

We first describe the pure carbon simulation, run cl 
in Table |T] We start by evolving the 27648 ion system 
at constant temperature for a time of a few million fm/c. 
During this time, we carefully adjust the temperature, by 
rescaling the velocities, so that about half of the system 
is solid and half is liquid. Then we evolve the system at 
constant energy for 50 million fm/c. Finally, as long as 
the potential is independent of momentum, the expec- 
tation value of the kinetic energy per particle is 3T/2. 
Therefore we estimate the melting temperature from the 



TABLE I: Computer Simulations with 27648 ions. The car- 
bon number fraction for the whole system is Xc, the total 
simulation time is t, and the final temperature is T. The car- 
bon number fraction of the solid phase is x^, while x^, is the 
carbon fraction of the liquid phase. 

Run x^ t (fm/c) T (MeV) x"^ x\ 

cl 1 5 X 10^ 0.02050(3) 1 1 

c90 0.90 8 X 10* 0.0195(1) 0.906(11) 0.900(8) 
c824 0.824 1 x 10^ 0.0197(1) 0.834(5) 0.819(5) 
c75 0.75 2 X 10^ 0.0193(1) 0.727(5) 0.766(5) 
c50 0.50 2 X 10^ 0.0217(1) 0.459(5) 0.525(5) 
c25 0.25 2 X 10^ 0.0256(1) 0.192(4) 0.284(4) 



kinetic energy. This yields a melting F value of 
F„ = 178.4 ±0.2, 



(2) 



see Table |l) This value is slightly larger than the F,„ = 
175 expected for the One Component Plasma (OCP) 
[TS][2D] and consistent with the value F^ « 178.2 pre- 
dicted by Eq. (4) in ref. [21] for a Yukawa system with 
our value of the parameter n = Xjin^l'^X) — 0.542 (and 
assuming F™ = 175 for the OCP). 

Our two-phase simulation method is subject to sys- 
tematic errors from finite size effects and from non- 
equilibrium effects due to the finite run time. Finite size 
effects could be important because the simulation vol- 
ume contains two liquid-solid interfaces so that an ion 
may be relatively close to one of the interfaces. Simula- 
tions in ref. [13] with 3456 ions yielded a melting F that 
is only slightly smaller F^ — 176.1 ± 0.7. We conclude 
that finite size effects are relatively small for our 27648 
ion single component system. 

We search for non-equilibrium effects by evaluating the 
temperature at times of 10, 20, 30, 40, and 50 million 
fm/c. We find very little time dependence. Therefore, we 
expect non-equilibrium effects to be small for our single 
component system. 

Our simulations for carbon / oxygen mixtures are per- 
formed in a similar way. We start from an initial con- 
figuration that is half liquid and half solid. The initial 
number fraction of carbon Xc in the solid phase is equal 
to that in the liquid phase. The system is evolved in 
time, first at constant temperature and then at constant 
energy. The carbon and oxygen ions are free to diffuse 
across the liquid-solid interfaces so that the number frac- 
tion of carbon in the liquid x\ can become different from 
the number fraction in the solid x%. 

We now present results for the phase diagram of carbon 
and oxygen mixtures. We measure the composition of the 
liquid x\ and solid x\ as follows. We divide the simulation 
volume into fifteen regions equally spaced in z coordinate 
and calculate the average (Z) for each region. Adjacent 
regions where (Z) changes from below to above average 
are assumed to represent liquid-solid interfaces and their 
composition is discarded. The remaining 11 regions are 
used to calculate the average composition of the liquid 
x\ and solid x*. Figure p] shows the difference x\ — x^ 



versus simulation time. Note that there are relatively 
large statistical errors. The solid curves in Fig. [l] are 
least squares fits of the form x^. — x^ = a[l — cxp{—bt)] 
with a and b constants. 
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FIG. 1: (Color on line) Number fraction of carbon in the 
liquid phase minus the number fraction of carbon in the solid 
phase versus simulation time for the simulations of Table [I] 



We see that non-equilibrium effects can be significantly 
larger for mixtures, because it can take a long time for 
impurities to diffuse. We note that the solid is enriched 
in oxygen, compared to the average composition, for runs 
c75, c50, and c25, while for run c824 the solid is enriched 
in carbon, compared to the average composition. Simple 
estimates of diffusion times suggest that the concentra- 
tion should have equilibrated by the relatively long sim- 
ulation time of 2 X 10^ fm/c. Runs c75, c50, and c25 
were performed on special purpose MDGRAPE-2 hard- 
ware J19j and took approximately six months of computer 
time each. 

We average the liquid and solid compositions over the 
final approximately 400 million fm/c of simulation time 
to determine the carbon-oxygen phase diagram. These 
results are listed in Table |l] and plotted in Fig. [2] as filled 
red circles. The upper curve gives the composition of the 
liquid that is in equilibrium with a solid of composition 
given by the lower curve. Note that the run cl is plot- 
ted twice in Fig. [2J first at Xq = (pure carbon) and 
then rescaled to a;o = 1 (pure oxygen). We find that the 
melting temperature of carbon oxygen mixtures is con- 
siderably below the constant F^ = 178 prediction, Eqs. 
|1|2[ This is plotted as a dot-dot-dashed line in Fig. [2] 
Our melting temperatures are also below the results of 
Segretain et al. [B]. We speculate that this could be 
because of small errors in Segretain et al's density func- 
tional calculations of the solid free energy. 

Our results agree qualitatively with Ogata et al. [7], 
and, in general, agree well with Medin and Gumming 
1131. Both of these calculations are based on Monte Carlo 
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FIG. 2; (Color on line) Melting temperature T of carbon / 
oxygen mixtures over the melting temperature Tc of pure car- 
bon, versus oxygen number fraction Xo = 1 — Xc- Simulation 
results of Table IT] are plotted as filled red circles connected 
by dashed lines. The Xo in the liquid and in the solid are 
shown as two separate lines. Also shown are the phase dia- 
gram results of Medin and Gumming [13] as solid black lines, 
the Ogata et al. results [Tj as dotted blue lines and the Seg- 
retain et al. results [6] as dot-dashed green lines. Finally the 
black dot-dot-dashed line corresponds to F = 178.4 in Eq. [l] 



or MD simulation free energies for the liquid and solid 
phases. Although the overall agreement with Medin and 
Gumming is good, there is a tendency for our simulations 
to predict smaller differences in composition between the 



liquid and solid phases x^ 



This could be because 



of finite size effects in our simulations. In equilibrium, 
there is a composition gradient, as a function of posi- 
tion, across the liquid-solid interface. Therefore, if one 
probes the composition of the liquid and solid in posi- 
tions that are too close to the interface, one will natu- 
rally get smaller differences between xj, and x'^. Alter- 



natively, Xq 



could be sensitive to small errors in 



Medin and Gummings' free energies. Overall, given the 
agreement between our results and those of Ogata et al. 
and Medin and Gumming, we conclude that the carbon 
/ oxygen phase diagram is largely known and that it is 
of azeotrope, instead of spindle, form. 

We now discuss implications of our carbon oxygen 
phase diagram on White Dwarf (WD) star crystalliza- 
tion, limits on the oxygen fraction of WDs, and possi- 
ble limits on the ^^G(a,7) reaction rate. Winget et al. 
observe the luminosity function (number of WD with a 
given luminosity versus luminosity) for the globular star 
cluster NGC 6397 [^ . They find a peak in the luminosity 
function that they attribute to crystallization, and they 
claim that the location of the peak is sensitive to the 
crystallization temperature of the WD core. 

Winget et al.'s observations agree well with theoreti- 



cal luminosity functions for 0.5-0.535 Mq WDs with pure 
carbon cores. The observations disagree with a theoreti- 
cal luminosity function assuming a WD core of 50% car- 
bon and 50% oxygen by mass (or Xo = 0.43 by number). 
This luminosity function fixed the melting temperature 
with Eqs. 
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^_^ 22J, for which T = 1.26Tc at Xo = 0.43, 
see the dot-dot-dashed line in Fig. [2] 

For simplicity, we assume theoretical luminosity func- 
tions can be characterized by the melting temperature of 
the core. The data clearly favor T near Tc and strongly 
disfavor T = 1.26Tc. If, in the future, one could set a 
limit of, for example, half this difference, T < 1.1 3Tc 
by comparing theoretical luminosity functions to obser- 
vations, then we can place limits on core oxygen concen- 
trations. These limits follow from the form of our phase 
diagram in Fig. [2J The melting temperature is close to Tc 
for Xo < 0.5 and then rises rapidly with increasing oxy- 
gen concentration. If one had a constraint of T < 1.13rc 
then our phase diagram implies 



<0.57 



(3) 



for the oxygen concentration by number or Xo < 0.64 
for the oxygen concentration by mass. We conclude that 
constraining the melting temperature of WD cores to be 
close to that for pure carbon constrains the oxygen con- 
centration to be of order 60% or less. 

Salaris et al. find that the oxygen concentration in WD 
cores depends on the cross section for the ^^C(a,7)^^0 
reaction, that can described by the astrophysical S fac- 
tor, and on the treatment of convection in stellar evolu- 
tion models [23j . With their treatment of convection and 
an effective 5* factor, at an energy of 300 keV, of ^300 — 
240 keV barns, Salaris et al. find Xo = 0.79 for the oxy- 
gen concentration in the core of a 0.54 Mq WD. This 
oxygen concentration would be ruled out if T < LIST^,. 



Alternatively, Salaris et al. find Xo = 0.57 for a 0.6 Mq 
WD if they assume a smaller value 5*300 = 170 keV barns. 
We expect the central oxygen abundance of a 0.54 Mq 
star to be slightly larger than that for a 0.6 Mq star and 
close to our limit of Xo < 0.64. We conclude that assum- 
ing T < 1.13rc, along with the Salaris et al. assumptions 
for convection, implies a limit on the effective S factor of 
order 



^300 < 170 keV barns. 



(4) 



This limit is consistent with the recent experimental de- 
termination by Buchmann and Barnes of 6*300 = 145 keV 
b, with an error of 25% to 35%, [1]. It is also consistent 
with the results of Tur, Heger and Austin who evaluate 
nucleosynthesis yields by varying both the triple-alpha 
and ^^C(a, 7) rates [21]. Their best fit value is 6*300 = 174 
keV b, with perhaps significant error. 

In the future, WD luminosity functions should be cal- 
culated using our carbon oxygen phase diagram. Possible 
limits on crystallization temperatures should be deduced 
from observations including careful consideration of sys- 
tematic errors. Finally, we will study the role of a small 
neon abundance on the phase diagram and WD crystal- 
lization using three component MD simulations, see for 
example |25j . 
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